home *** CD-ROM | disk | FTP | other *** search
/ Enter 2004 January / enter-2004-01.iso / files / maxima-5.9.0.exe / {app} / share / maxima / 5.9.0 / src / numerical / slatec / zasyi.lisp < prev    next >
Encoding:
Text File  |  2003-02-09  |  9.1 KB  |  245 lines
  1. ;;; Compiled by f2cl version 2.0 beta 2002-05-06
  2. ;;; 
  3. ;;; Options: ((:prune-labels nil) (:auto-save t) (:relaxed-array-decls t)
  4. ;;;           (:coerce-assigns :as-needed) (:array-type ':simple-array)
  5. ;;;           (:array-slicing nil) (:declare-common nil)
  6. ;;;           (:float-format double-float))
  7.  
  8. (in-package "SLATEC")
  9.  
  10.  
  11. (let ((pi_ 3.141592653589793)
  12.       (rtpi 0.15915494309189535)
  13.       (zeror 0.0)
  14.       (zeroi 0.0)
  15.       (coner 1.0)
  16.       (conei 0.0))
  17.   (declare (type double-float conei coner zeroi zeror rtpi pi_))
  18.   (defun zasyi (zr zi fnu kode n yr yi nz rl tol elim alim)
  19.     (declare (type (simple-array double-float (*)) yr yi)
  20.              (type f2cl-lib:integer4 kode n nz)
  21.              (type double-float zr zi fnu rl tol elim alim))
  22.     (prog ((i 0) (ib 0) (il 0) (inu 0) (j 0) (jl 0) (k 0) (koded 0) (m 0)
  23.            (nn 0) (aa 0.0) (aez 0.0) (ak 0.0) (ak1i 0.0) (ak1r 0.0) (arg 0.0)
  24.            (arm 0.0) (atol 0.0) (az 0.0) (bb 0.0) (bk 0.0) (cki 0.0) (ckr 0.0)
  25.            (cs1i 0.0) (cs1r 0.0) (cs2i 0.0) (cs2r 0.0) (czi 0.0) (czr 0.0)
  26.            (dfnu 0.0) (dki 0.0) (dkr 0.0) (dnu2 0.0) (ezi 0.0) (ezr 0.0)
  27.            (fdn 0.0) (p1i 0.0) (p1r 0.0) (raz 0.0) (rtr1 0.0) (rzi 0.0)
  28.            (rzr 0.0) (s 0.0) (sgn 0.0) (sqk 0.0) (sti 0.0) (str 0.0) (s2i 0.0)
  29.            (s2r 0.0) (tzi 0.0) (tzr 0.0))
  30.       (declare
  31.        (type double-float tzr tzi s2r s2i str sti sqk sgn s rzr rzi rtr1 raz
  32.         p1r p1i fdn ezr ezi dnu2 dkr dki dfnu czr czi cs2r cs2i cs1r cs1i ckr
  33.         cki bk bb az atol arm arg ak1r ak1i ak aez aa)
  34.        (type f2cl-lib:integer4 nn m koded k jl j inu il ib i))
  35.       (setf nz 0)
  36.       (setf az (zabs zr zi))
  37.       (setf arm (* 1000.0 (f2cl-lib:d1mach 1)))
  38.       (setf rtr1 (f2cl-lib:fsqrt arm))
  39.       (setf il (min (the f2cl-lib:integer4 2) (the f2cl-lib:integer4 n)))
  40.       (setf dfnu (+ fnu (f2cl-lib:int-sub n il)))
  41.       (setf raz (/ 1.0 az))
  42.       (setf str (* zr raz))
  43.       (setf sti (* (- zi) raz))
  44.       (setf ak1r (* rtpi str raz))
  45.       (setf ak1i (* rtpi sti raz))
  46.       (multiple-value-bind
  47.           (var-0 var-1 var-2 var-3)
  48.           (zsqrt ak1r ak1i ak1r ak1i)
  49.         (declare (ignore var-0 var-1))
  50.         (setf ak1r var-2)
  51.         (setf ak1i var-3))
  52.       (setf czr zr)
  53.       (setf czi zi)
  54.       (if (/= kode 2) (go label10))
  55.       (setf czr zeror)
  56.       (setf czi zi)
  57.      label10
  58.       (if (> (abs czr) elim) (go label100))
  59.       (setf dnu2 (+ dfnu dfnu))
  60.       (setf koded 1)
  61.       (if (and (> (abs czr) alim) (> n 2)) (go label20))
  62.       (setf koded 0)
  63.       (multiple-value-bind
  64.           (var-0 var-1 var-2 var-3)
  65.           (zexp czr czi str sti)
  66.         (declare (ignore var-0 var-1))
  67.         (setf str var-2)
  68.         (setf sti var-3))
  69.       (multiple-value-bind
  70.           (var-0 var-1 var-2 var-3 var-4 var-5)
  71.           (zmlt ak1r ak1i str sti ak1r ak1i)
  72.         (declare (ignore var-0 var-1 var-2 var-3))
  73.         (setf ak1r var-4)
  74.         (setf ak1i var-5))
  75.      label20
  76.       (setf fdn 0.0)
  77.       (if (> dnu2 rtr1) (setf fdn (* dnu2 dnu2)))
  78.       (setf ezr (* zr 8.0))
  79.       (setf ezi (* zi 8.0))
  80.       (setf aez (* 8.0 az))
  81.       (setf s (/ tol aez))
  82.       (setf jl (f2cl-lib:int (+ rl rl 2)))
  83.       (setf p1r zeror)
  84.       (setf p1i zeroi)
  85.       (if (= zi 0.0) (go label30))
  86.       (setf inu (f2cl-lib:int fnu))
  87.       (setf arg (* (- fnu inu) pi_))
  88.       (setf inu (f2cl-lib:int-sub (f2cl-lib:int-add inu n) il))
  89.       (setf ak (- (sin arg)))
  90.       (setf bk (cos arg))
  91.       (if (< zi 0.0) (setf bk (- bk)))
  92.       (setf p1r ak)
  93.       (setf p1i bk)
  94.       (if (= (mod inu 2) 0) (go label30))
  95.       (setf p1r (- p1r))
  96.       (setf p1i (- p1i))
  97.      label30
  98.       (f2cl-lib:fdo (k 1 (f2cl-lib:int-add k 1))
  99.                     ((> k il) nil)
  100.         (tagbody
  101.           (setf sqk (- fdn 1.0))
  102.           (setf atol (* s (abs sqk)))
  103.           (setf sgn 1.0)
  104.           (setf cs1r coner)
  105.           (setf cs1i conei)
  106.           (setf cs2r coner)
  107.           (setf cs2i conei)
  108.           (setf ckr coner)
  109.           (setf cki conei)
  110.           (setf ak 0.0)
  111.           (setf aa 1.0)
  112.           (setf bb aez)
  113.           (setf dkr ezr)
  114.           (setf dki ezi)
  115.           (f2cl-lib:fdo (j 1 (f2cl-lib:int-add j 1))
  116.                         ((> j jl) nil)
  117.             (tagbody
  118.               (multiple-value-bind
  119.                   (var-0 var-1 var-2 var-3 var-4 var-5)
  120.                   (zdiv ckr cki dkr dki str sti)
  121.                 (declare (ignore var-0 var-1 var-2 var-3))
  122.                 (setf str var-4)
  123.                 (setf sti var-5))
  124.               (setf ckr (* str sqk))
  125.               (setf cki (* sti sqk))
  126.               (setf cs2r (+ cs2r ckr))
  127.               (setf cs2i (+ cs2i cki))
  128.               (setf sgn (- sgn))
  129.               (setf cs1r (+ cs1r (* ckr sgn)))
  130.               (setf cs1i (+ cs1i (* cki sgn)))
  131.               (setf dkr (+ dkr ezr))
  132.               (setf dki (+ dki ezi))
  133.               (setf aa (/ (* aa (abs sqk)) bb))
  134.               (setf bb (+ bb aez))
  135.               (setf ak (+ ak 8.0))
  136.               (setf sqk (- sqk ak))
  137.               (if (<= aa atol) (go label50))
  138.              label40))
  139.           (go label110)
  140.          label50
  141.           (setf s2r cs1r)
  142.           (setf s2i cs1i)
  143.           (if (>= (+ zr zr) elim) (go label60))
  144.           (setf tzr (+ zr zr))
  145.           (setf tzi (+ zi zi))
  146.           (multiple-value-bind
  147.               (var-0 var-1 var-2 var-3)
  148.               (zexp (- tzr) (- tzi) str sti)
  149.             (declare (ignore var-0 var-1))
  150.             (setf str var-2)
  151.             (setf sti var-3))
  152.           (multiple-value-bind
  153.               (var-0 var-1 var-2 var-3 var-4 var-5)
  154.               (zmlt str sti p1r p1i str sti)
  155.             (declare (ignore var-0 var-1 var-2 var-3))
  156.             (setf str var-4)
  157.             (setf sti var-5))
  158.           (multiple-value-bind
  159.               (var-0 var-1 var-2 var-3 var-4 var-5)
  160.               (zmlt str sti cs2r cs2i str sti)
  161.             (declare (ignore var-0 var-1 var-2 var-3))
  162.             (setf str var-4)
  163.             (setf sti var-5))
  164.           (setf s2r (+ s2r str))
  165.           (setf s2i (+ s2i sti))
  166.          label60
  167.           (setf fdn (+ fdn (* 8.0 dfnu) 4.0))
  168.           (setf p1r (- p1r))
  169.           (setf p1i (- p1i))
  170.           (setf m (f2cl-lib:int-add (f2cl-lib:int-sub n il) k))
  171.           (f2cl-lib:fset (f2cl-lib:fref yr (m) ((1 n)))
  172.                          (- (* s2r ak1r) (* s2i ak1i)))
  173.           (f2cl-lib:fset (f2cl-lib:fref yi (m) ((1 n)))
  174.                          (+ (* s2r ak1i) (* s2i ak1r)))
  175.          label70))
  176.       (if (<= n 2) (go end_label))
  177.       (setf nn n)
  178.       (setf k (f2cl-lib:int-sub nn 2))
  179.       (setf ak (coerce (the f2cl-lib:integer4 k) 'double-float))
  180.       (setf str (* zr raz))
  181.       (setf sti (* (- zi) raz))
  182.       (setf rzr (* (+ str str) raz))
  183.       (setf rzi (* (+ sti sti) raz))
  184.       (setf ib 3)
  185.       (f2cl-lib:fdo (i ib (f2cl-lib:int-add i 1))
  186.                     ((> i nn) nil)
  187.         (tagbody
  188.           (f2cl-lib:fset (f2cl-lib:fref yr (k) ((1 n)))
  189.                          (+
  190.                           (* (+ ak fnu)
  191.                              (-
  192.                               (* rzr
  193.                                  (f2cl-lib:fref yr
  194.                                                 ((f2cl-lib:int-add k 1))
  195.                                                 ((1 n))))
  196.                               (* rzi
  197.                                  (f2cl-lib:fref yi
  198.                                                 ((f2cl-lib:int-add k 1))
  199.                                                 ((1 n))))))
  200.                           (f2cl-lib:fref yr ((f2cl-lib:int-add k 2)) ((1 n)))))
  201.           (f2cl-lib:fset (f2cl-lib:fref yi (k) ((1 n)))
  202.                          (+
  203.                           (* (+ ak fnu)
  204.                              (+
  205.                               (* rzr
  206.                                  (f2cl-lib:fref yi
  207.                                                 ((f2cl-lib:int-add k 1))
  208.                                                 ((1 n))))
  209.                               (* rzi
  210.                                  (f2cl-lib:fref yr
  211.                                                 ((f2cl-lib:int-add k 1))
  212.                                                 ((1 n))))))
  213.                           (f2cl-lib:fref yi ((f2cl-lib:int-add k 2)) ((1 n)))))
  214.           (setf ak (- ak 1.0))
  215.           (setf k (f2cl-lib:int-sub k 1))
  216.          label80))
  217.       (if (= koded 0) (go end_label))
  218.       (multiple-value-bind
  219.           (var-0 var-1 var-2 var-3)
  220.           (zexp czr czi ckr cki)
  221.         (declare (ignore var-0 var-1))
  222.         (setf ckr var-2)
  223.         (setf cki var-3))
  224.       (f2cl-lib:fdo (i 1 (f2cl-lib:int-add i 1))
  225.                     ((> i nn) nil)
  226.         (tagbody
  227.           (setf str
  228.                   (- (* (f2cl-lib:fref yr (i) ((1 n))) ckr)
  229.                      (* (f2cl-lib:fref yi (i) ((1 n))) cki)))
  230.           (f2cl-lib:fset (f2cl-lib:fref yi (i) ((1 n)))
  231.                          (+ (* (f2cl-lib:fref yr (i) ((1 n))) cki)
  232.                             (* (f2cl-lib:fref yi (i) ((1 n))) ckr)))
  233.           (f2cl-lib:fset (f2cl-lib:fref yr (i) ((1 n))) str)
  234.          label90))
  235.       (go end_label)
  236.      label100
  237.       (setf nz -1)
  238.       (go end_label)
  239.      label110
  240.       (setf nz -2)
  241.       (go end_label)
  242.      end_label
  243.       (return (values nil nil nil nil nil nil nil nz nil nil nil nil)))))
  244.  
  245.